############### R-CODE1: RUN ONLY AFTER RUNNING PART1 of STATA CODE ###############


###Figure 6: First-Period and Last-Period Contributions by City (Herrmann et al., 2008)

### constructing herrmann 1st-10th period data
#install plyr package
#install.packages("plyr") 
library(plyr)
#install.packages("ggplot2")
library(ggplot2)
# set working directory to corresponding main folder
setwd("/Users/mustafayahsi/Desktop/The Seeds of Success")
allcitiesp <- read.csv("Stata/Stata Produced Data/herrmannallcitiesp.csv")
a<-allcitiesp
### Rearranging data so that 1st & 10th period avg cont are on the same row
a$city<-as.factor(a$city)
numrow<-nrow(a)/2
b<-array(NA,dim=c(numrow,5))
j<-1
k<-count(a$city=="Athens")[2,2]/2
for(i in j:k){
  b[i,1]=a$secseq[2*(i-1)+1]
  b[i,2]=a$group_avg[2*(i-1)+1]
  b[i,3]=a$group_avg[2*i]
  b[i,5]=a$grp[2*(i-1)+1]
  b[i,4]=1
}
j<-k+1
k<-k+count(a$city=="Bonn")[2,2]/2
for(i in j:k){
  b[i,1]=a$secseq[2*(i-1)+1]
  b[i,2]=a$group_avg[2*(i-1)+1]
  b[i,3]=a$group_avg[2*i]
  b[i,5]=a$grp[2*(i-1)+1]
  b[i,4]=2
}

j<-k+1
k<-k+count(a$city=="Boston")[2,2]/2
for(i in j:k){
  b[i,1]=a$secseq[2*(i-1)+1]
  b[i,2]=a$group_avg[2*(i-1)+1]
  b[i,3]=a$group_avg[2*i]
  b[i,5]=a$grp[2*(i-1)+1]
  b[i,4]=3
  
}

j<-k+1
k<-k+count(a$city=="Chengdu")[2,2]/2
for(i in j:k){
  b[i,1]=a$secseq[2*(i-1)+1]
  b[i,2]=a$group_avg[2*(i-1)+1]
  b[i,3]=a$group_avg[2*i]
  b[i,5]=a$grp[2*(i-1)+1]
  b[i,4]=4
  
}
j<-k+1
k<-k+count(a$city=="Copenhagen")[2,2]/2
for(i in j:k){
  b[i,1]=a$secseq[2*(i-1)+1]
  b[i,2]=a$group_avg[2*(i-1)+1]
  b[i,3]=a$group_avg[2*i]
  b[i,5]=a$grp[2*(i-1)+1]
  b[i,4]=5
  
}
j<-k+1
k<-k+count(a$city=="Dnipropetrovs'k")[2,2]/2
for(i in j:k){
  b[i,1]=a$secseq[2*(i-1)+1]
  b[i,2]=a$group_avg[2*(i-1)+1]
  b[i,3]=a$group_avg[2*i]
  b[i,5]=a$grp[2*(i-1)+1]
  b[i,4]=6
  
}
j<-k+1
k<-k+count(a$city=="Istanbul")[2,2]/2
for(i in j:k){
  b[i,1]=a$secseq[2*(i-1)+1]
  b[i,2]=a$group_avg[2*(i-1)+1]
  b[i,3]=a$group_avg[2*i]
  b[i,5]=a$grp[2*(i-1)+1]
  b[i,4]=7
  
}
j<-k+1
k<-k+count(a$city=="Melbourne")[2,2]/2
for(i in j:k){
  b[i,1]=a$secseq[2*(i-1)+1]
  b[i,2]=a$group_avg[2*(i-1)+1]
  b[i,3]=a$group_avg[2*i]
  b[i,5]=a$grp[2*(i-1)+1]
  b[i,4]=8
  
}
j<-k+1
k<-k+count(a$city=="Minsk")[2,2]/2
for(i in j:k){
  b[i,1]=a$secseq[2*(i-1)+1]
  b[i,2]=a$group_avg[2*(i-1)+1]
  b[i,3]=a$group_avg[2*i]
  b[i,5]=a$grp[2*(i-1)+1]
  b[i,4]=9
  
}
j<-k+1
k<-k+count(a$city=="Muscat")[2,2]/2
for(i in j:k){
  b[i,1]=a$secseq[2*(i-1)+1]
  b[i,2]=a$group_avg[2*(i-1)+1]
  b[i,3]=a$group_avg[2*i]
  b[i,5]=a$grp[2*(i-1)+1]
  b[i,4]=10
  
}
j<-k+1
k<-k+count(a$city=="Nottingham")[2,2]/2
for(i in j:k){
  b[i,1]=a$secseq[2*(i-1)+1]
  b[i,2]=a$group_avg[2*(i-1)+1]
  b[i,3]=a$group_avg[2*i]
  b[i,5]=a$grp[2*(i-1)+1]
  b[i,4]=11
  
}
j<-k+1
k<-k+count(a$city=="Riyadh")[2,2]/2
for(i in j:k){
  b[i,1]=a$secseq[2*(i-1)+1]
  b[i,2]=a$group_avg[2*(i-1)+1]
  b[i,3]=a$group_avg[2*i]
  b[i,5]=a$grp[2*(i-1)+1]
  b[i,4]=12
  
}
j<-k+1
k<-k+count(a$city=="Samara")[2,2]/2
for(i in j:k){
  b[i,1]=a$secseq[2*(i-1)+1]
  b[i,2]=a$group_avg[2*(i-1)+1]
  b[i,3]=a$group_avg[2*i]
  b[i,5]=a$grp[2*(i-1)+1]
  b[i,4]=13
  
}
j<-k+1
k<-k+count(a$city=="Seoul")[2,2]/2
for(i in j:k){
  b[i,1]=a$secseq[2*(i-1)+1]
  b[i,2]=a$group_avg[2*(i-1)+1]
  b[i,3]=a$group_avg[2*i]
  b[i,5]=a$grp[2*(i-1)+1]
  b[i,4]=14
  
}
j<-k+1
k<-k+count(a$city=="St. Gallen")[2,2]/2
for(i in j:k){
  b[i,1]=a$secseq[2*(i-1)+1]
  b[i,2]=a$group_avg[2*(i-1)+1]
  b[i,3]=a$group_avg[2*i]
  b[i,5]=a$grp[2*(i-1)+1]
  b[i,4]=15
  
}
j<-k+1
k<-k+count(a$city=="Zurich")[2,2]/2
for(i in j:k){
  b[i,1]=a$secseq[2*(i-1)+1]
  b[i,2]=a$group_avg[2*(i-1)+1]
  b[i,3]=a$group_avg[2*i]
  b[i,5]=a$grp[2*(i-1)+1]
  b[i,4]=16
  
}

## Rename columns
colnames(b)<-c("secseq","group_avg1","group_avg2","City","group")
b<-as.data.frame(b)
b$City<-as.factor(b$City)
b$secseq<-as.factor(b$secseq)
b$city2<-b$City
b$City<-revalue(b$City, c("1"="Athens", "2"="Bonn","3"="Boston", "4"="Chengdu","5"="Copenhagen", "6"="Dnipropetrovs'k","7"="Istanbul",
                          "8"="Melbourne","9"="Minsk", "10"="Muscat","11"="Nottingham", "12"="Riyadh","13"="Samara", 
                          "14"="Seoul","15"="St. Gallen", "16"="Zurich"))

#### Figure 6: First-Period and Last-Period Contributions by City (Herrmann et al., 2008) ######
b2<-subset(b,b$City=="Minsk" | b$City=="Muscat" | b$City=="Riyadh" | b$City=="Samara" | b$City=="Athens" | b$City=="Dnipropetrovs'k")
b2<-as.data.frame(b2)
b2$City<-as.factor(b2$City)
p1<-ggplot(b2,aes(x=group_avg1,y=group_avg2))+facet_wrap(~City)+labs(y = "Contribution - 10th Period",x="Contribution - 1st Period")+theme(
  strip.background = element_blank(),
  strip.text.x = element_blank()
)+scale_y_continuous(breaks = c(0,5,10,15,20),label = c("0","5","10","15","20"))+
  scale_x_continuous(limits=c(0, 20),breaks = c(0,5,10,15,20),label = c("0","5","10","15","20"))+scale_shape_manual(values=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16))
p2<-p1+geom_point()+theme(strip.background = element_blank(), strip.text.x = element_blank())+scale_color_discrete(guide=FALSE)
p3<-p2+geom_abline(slope=1,intercept=0,linetype="dashed")+geom_smooth(method='lm',formula=y~x,se=TRUE,color="black")
p4<-p3+theme(strip.background = element_blank(), strip.text.x = element_blank())
p5<-p4+theme_bw()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p5+ theme(legend.position="bottom")


## write to csv file to prepare for PART2 CODE OF STATA*******
write.csv(b,"R-Studio/R Produced Data/allcities.csv")
write.csv(b2,"R-Studio/R Produced Data/6cities.csv")

############################# R CODE2: RUN ONLY AFTER RUNNING PART 2 STATA CODE #############
library(plyr)
library(ggplot2)
sixcitiesp <- read.csv("Stata/Stata Produced Data/herrmann-fe6cities.csv")
d<-sixcitiesp

##### Figure D.2: First-Period and Last-Period Contributions by City (FE) (Herrmann et al., 2008) ######

p1<-ggplot(d,aes(x=group_avg1,y=group_avg2,shape=city))+labs(y = "Contribution - 10th Period",x="Contribution - 1st Period")+theme(
  strip.background = element_blank(),
  strip.text.x = element_blank()
)+scale_y_continuous(breaks = c(0,5,10,15,20),label = c("0","5","10","15","20"))+
  scale_x_continuous(limits=c(0, 20),breaks = c(0,5,10,15,20),label = c("0","5","10","15","20"))+scale_shape_manual(values=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16))
p2<-p1+geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci),fill="grey")
p3<-p2+geom_point()+theme(strip.background = element_blank(), strip.text.x = element_blank())+scale_color_discrete(guide=FALSE)
p4<-p3+geom_abline(slope=1,intercept=0,linetype="dashed")
p5<-p4+geom_line(aes(y=fitted))
p6<-p5+theme(strip.background = element_blank(), strip.text.x = element_blank())
p7<-p6+theme_bw()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p7+ theme(legend.position="bottom")



### Figure 5: Comparison of Average Contributions in the P-experiment ######
#### herrmann part of the graph
b3<-subset(b,b$City=="Istanbul" )
b3<-as.data.frame(b3)
b3$City<-as.factor(b3$City)
p1<-ggplot(b3,aes(x=group_avg1,y=group_avg2))+labs(y = "Contribution - 10th Period",x="Contribution - 1st Period")+theme(
  strip.background = element_blank(),
  strip.text.x = element_blank()
)+scale_y_continuous(breaks = c(0,5,10,15,20),label = c("0","5","10","15","20"))+
  scale_x_continuous(limits=c(0, 20),breaks = c(0,5,10,15,20),label = c("0","5","10","15","20"))+scale_shape_manual(values=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16))
p2<-p1+geom_point()+theme(strip.background = element_blank(), strip.text.x = element_blank())+scale_color_discrete(guide=FALSE)
p3<-p2+geom_abline(slope=1,intercept=0,linetype="dashed")+geom_smooth(method='lm',formula=y~x,se=TRUE,color="black")
p4<-p3+theme(strip.background = element_blank(), strip.text.x = element_blank())
p5<-p4+theme_bw()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p5+ theme(legend.position="bottom")

########### Figure 5: Our Data Part of the Graph 
#### Constructing trueavg for our data 
#install.packages(readxl)
#install.packages(openxlsx)
library(readxl)
library(openxlsx)
k<-read_excel("Stata/Stata Produced Data/oursaverage.xlsx",col_names = FALSE)
k<-as.data.frame(k)
tavg<-matrix(NA,15,10)
for(i in 1:10){
  for(j in 1:15){
    tavg[j,i]=k[j+15*(i-1),3]
  }
}
tavg_df <- as.data.frame(tavg)
write.xlsx(tavg_df,file="R-Studio/R produced data/trueavg.xlsx")

### read true average contributions from data
library(readxl)
trueavg<- read_excel("R-Studio/R produced data/trueavg.xlsx") 
trueavg<-as.matrix(trueavg)
b4<-matrix(NA,15,2)
b4[,1]<-trueavg[,1]
b4[,2]<-trueavg[,10]
b4<-as.data.frame(b4)
## Regress last period contribution on 1st period contribution
colnames(b4)<-c("group_avg1","group_avg2")
summary(lm(group_avg2 ~ group_avg1,data=b4))

#### Figure 5, second part
p1<-ggplot(b4,aes(x=group_avg1,y=group_avg2))+labs(y = "Contribution - 10th Period",x="Contribution - 1st Period")+theme(
  strip.background = element_blank(),
  strip.text.x = element_blank()
)+scale_y_continuous(breaks = c(0,5,10,15,20),label = c("0","5","10","15","20"))+
  scale_x_continuous(limits=c(0, 20),breaks = c(0,5,10,15,20),label = c("0","5","10","15","20"))+scale_shape_manual(values=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16))
#p2<-p1+geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci),fill="grey")
p3<-p1+geom_point()+theme(strip.background = element_blank(), strip.text.x = element_blank())+scale_color_discrete(guide=FALSE)
p4<-p3+geom_abline(slope=1,intercept=0,linetype="dashed")+geom_smooth(method='lm',formula=y~x,se=TRUE,color="black")
#p5<-p4+geom_line(aes(y=fitted))
p6<-p4+theme(strip.background = element_blank(), strip.text.x = element_blank())
p7<-p6+theme_bw()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
p7+ theme(legend.position="bottom")



########################### AGENT BASED SIMULATION PART########################



library(extrafont)
loadfonts(device="win")       #Register fonts for Windows bitmap output
par(family="Times")
##########################readexcel
library(readxl, quietly=TRUE)
firstperiod<- read_excel("Stata/Stata Produced Data/firstperiod.xlsx")
trueavg<- read_excel("R-Studio/R Produced Data/trueavg.xlsx")
trueavg<-as.matrix(trueavg)
head(firstperiod)
firstperiod<-as.data.frame(firstperiod)
data<-matrix(NA,15,4)
for(i in 1:15){
  for(j in 1:4){
    data[i,j]<-firstperiod$contribution[4*(i-1)+j]
  }
}
##############################betas
r=100
contribution<-array(NA,dim=c(r,15,4,21))
for(i in 1:15){
  for (z in 1:100){
    contribution[z,i,,1]=data[i,]
  }
}

probability<-array(NA,dim=c(r,50,4,4,21))
probabilitylogit<-array(NA,dim=c(r,50,4,21))
punishment<-array(NA,dim=c(r,50,4,4,21))
period<-array(1:20)
contributionlogit<-array(NA,dim=c(r,50,4,21))
avgcontofthers<-array(NA,dim=c(r,50,4,21))
avgcontofremaining<-array(NA,dim=c(r,50,4,4,21))
avgcont<-array(NA,dim=c(r,50,21))
receivedpoints<-array(NA,dim=c(r,50,4,21))
sanction<-array(NA,dim=c(r,50,4,4,21))
increaselogit<-array(NA,dim=c(r,50,4,21))
increase<-array(NA,dim=c(r,50,4,21))
decreaselogit<-array(NA,dim=c(r,50,4,21))
decrease<-array(NA,dim=c(r,50,4,21))
averagecontribution<-array(NA,dim=c(r,15,20))
##########################definitions

############error term for contribution clustered groups
## retrieved from stata regressions 
sigmainc<- 0.3023655 
sigmadec<-0.2872764    

mean<-0
x1<-matrix(NA,4,10)
epsiloninc<-rnorm(x1,mean,sigmainc)
epsilondec<-rnorm(x1,mean,sigmadec)
epsiloninc<-array(epsiloninc,dim=c(r,50,4,21))
epsilondec<-array(epsilondec,dim=c(r,50,4,21))




######################################increase&decrease logit
## all coefficients are retrieved from relevant tables in the paper 
#### with no cons
ilreceived=0.175
dlreceived=0
dlothersavg=-0.318
dlcontprev=0.216
###ols reg cont
increcpoints=0.399
incavgcontofthers=0.473
inccontprev=0.511
incperiod=0.151

decrecpoints=-0.513
decavgcontofthers=0.495
deccontprev=0.412

########
lrcont<--0.146
lscont<-0.033
lavgcontsanc<-0.094
lrpoints<-0.045
lperiod<--0.082
period<-(1:10)
period<-(1:20)
## use only first period contribution to forecast next period contributions for 20 consecutive periods
#######inc/dec main part
for(z in 1:100){
  crvseed <-.Random.seed[z]
  for (t in 1:15){
    for(i in 1:20){
      avgcont[z,t,i]=(sum(contribution[z,t,,i]))/4
      if(i==1){
        for(j in 1:4){
          avgcontofthers[z,t,j,i]<-(sum(contribution[z,t,,i])-contribution[z,t,j,i])/3
          for(k in 1:4){
            avgcontofremaining[z,t,j,k,i]<-(sum(contribution[z,t,,i])-contribution[z,t,j,i]-contribution[z,t,k,i])/2
            if(contribution[z,t,j,i]>contribution[z,t,k,i]){
              punishment[z,t,j,k,i]<--0.100*contribution[z,t,k,i]  + 0.067*avgcontofremaining[z,t,j,k,i]
              probability[z,t,j,k,i]<-1/(1+exp(-punishment[z,t,j,k,i]))
              sanction[z,t,j,k,i]<-rbinom(1,1,probability[z,t,j,k,i])
              if( k==j){
                sanction[z,t,j,k,i]=0
              }
            }
            else if(contribution[z,t,k,i]>=contribution[z,t,j,i]){
              punishment[z,t,j,k,i]<--0.085*contribution[z,t,k,i] 
              probability[z,t,j,k,i]<-1/(1+exp(-punishment[z,t,j,k,i]))
              sanction[z,t,j,k,i]<-rbinom(1,1,probability[z,t,j,k,i])
              if( k==j){
                sanction[z,t,j,k,i]=0
              }
            }
          }
        }
      }
      if(i!=1){
        for(j in 1:4){
          avgcontofthers[z,t,j,i]<-(sum(contribution[z,t,,i])-contribution[z,t,j,i])/3
          receivedpoints[z,t,j,i-1]<-sum(sanction[z,t,,j,i-1])
          for(k in 1:4){
            avgcontofremaining[z,t,j,k,i]<-(sum(contribution[z,t,,i])-contribution[z,t,j,i]-contribution[z,t,k,i])/2
            punishment[z,t,j,k,i]<-lrcont*contribution[z,t,k,i] + lscont*contribution[z,t,j,i] + lavgcontsanc*avgcontofremaining[z,t,j,k,i]+lrpoints*receivedpoints[z,t,j,i-1]+ lperiod*period[i]
            probability[z,t,j,k,i]<-1/(1+exp(-punishment[z,t,j,k,i]))
            sanction[z,t,j,k,i]<-rbinom(1,1,probability[z,t,j,k,i])
            if( k==j){
              sanction[z,t,j,k,i]=0
            }
          }
        }
      }
      for(m in 1:4){
        if(contribution[z,t,m,i]<avgcont[z,t,i]){
          avgcontofthers[z,t,m,i]<-(sum(contribution[z,t,,i])-contribution[z,t,m,i])/3
          receivedpoints[z,t,m,i]=sum(sanction[z,t,,m,i])
          increaselogit[z,t,m,i+1]=ilreceived*receivedpoints[z,t,m,i]
          probabilitylogit[z,t,m,i+1]<-1/(1+exp(-increaselogit[z,t,m,i+1]))
          increase[z,t,m,i+1]<-rbinom(1,1,probabilitylogit[z,t,m,i+1])
          if(increase[z,t,m,i+1]== 1){
            contribution[z,t,m,i+1]<-increcpoints*receivedpoints[z,t,m,i]+inccontprev*contribution[z,t,m,i]+incavgcontofthers*avgcontofthers[z,t,m,i]+incperiod*period[i]+epsiloninc[z,t,m,i+1]
          }
          else if(increase[z,t,m,i+1]== 0){
            contribution[z,t,m,i+1]<-contribution[z,t,m,i]
          }
          
        }
        else if (contribution[z,t,m,i]>=avgcont[z,t,i]){
          avgcontofthers[z,t,m,i]<-(sum(contribution[z,t,,i])-contribution[z,t,m,i])/3
          receivedpoints[z,t,m,i]=sum(sanction[z,t,,m,i])
          decreaselogit[z,t,m,i+1]=dlreceived*receivedpoints[z,t,m,i]+dlothersavg*avgcontofthers[z,t,m,i]+dlcontprev*contribution[z,t,m,i] ###+dlconstant
          probabilitylogit[z,t,m,i+1]<-1/(1+exp(-decreaselogit[z,t,m,i+1]))
          decrease[z,t,m,i+1]<-rbinom(1,1,probabilitylogit[z,t,m,i+1])
          if(decrease[z,t,m,i+1]== 1){
            contribution[z,t,m,i+1]<-decrecpoints*receivedpoints[z,t,m,i]+deccontprev*contribution[z,t,m,i]+decavgcontofthers*avgcontofthers[z,t,m,i]+epsilondec[z,t,m,i+1]
          }
          else if(decrease[z,t,m,i+1]== 0){
            contribution[z,t,m,i+1]<-contribution[z,t,m,i]
          }
          
        }
      }
    }
  }
  
  for (t in 1:15){
    for(i in 1:4){
      receivedpoints[z,t,i,20]=sum(sanction[z,t,,i,20])
    }
  }
  
  for (i in 1:15){
    for(j in 1:20){
      averagecontribution[z,i,j]=sum(contribution[z,i,,j])/4
      
    }
  }
}
## arrange forecasted data to construct relevant graphs and retrieve confidence intervals
averagecontribution.t<-array(NA,dim=c(r,20,15))
for(z in 1:100){
  averagecontribution.t[z,,]<-t(averagecontribution[z,,])
  
}
averagegroup1<-array(NA,dim=c(20,100))
averagegroup2<-array(NA,dim=c(20,100))
averagegroup3<-array(NA,dim=c(20,100))
averagegroup4<-array(NA,dim=c(20,100))
averagegroup5<-array(NA,dim=c(20,100))
averagegroup6<-array(NA,dim=c(20,100))
averagegroup7<-array(NA,dim=c(20,100))
averagegroup8<-array(NA,dim=c(20,100))
averagegroup9<-array(NA,dim=c(20,100))
averagegroup10<-array(NA,dim=c(20,100))
averagegroup11<-array(NA,dim=c(20,100))
averagegroup12<-array(NA,dim=c(20,100))
averagegroup13<-array(NA,dim=c(20,100))
averagegroup14<-array(NA,dim=c(20,100))
averagegroup15<-array(NA,dim=c(20,100))


for(i in 1:100){
  averagegroup1[,i]=averagecontribution.t[i,,1]
  averagegroup2[,i]=averagecontribution.t[i,,2]
  averagegroup3[,i]=averagecontribution.t[i,,3]
  averagegroup4[,i]=averagecontribution.t[i,,4]
  averagegroup5[,i]=averagecontribution.t[i,,5]
  averagegroup6[,i]=averagecontribution.t[i,,6]
  averagegroup7[,i]=averagecontribution.t[i,,7]
  averagegroup8[,i]=averagecontribution.t[i,,8]
  averagegroup9[,i]=averagecontribution.t[i,,9]
  averagegroup10[,i]=averagecontribution.t[i,,10]
  averagegroup11[,i]=averagecontribution.t[i,,11]
  averagegroup12[,i]=averagecontribution.t[i,,12]
  averagegroup13[,i]=averagecontribution.t[i,,13]
  averagegroup14[,i]=averagecontribution.t[i,,14]
  averagegroup15[,i]=averagecontribution.t[i,,15]
}


avg1=array(NA,dim=c(10,1))
avg2=array(NA,dim=c(10,1))
avg3=array(NA,dim=c(10,1))
avg4=array(NA,dim=c(10,1))
avg5=array(NA,dim=c(10,1))
avg6=array(NA,dim=c(10,1))
avg7=array(NA,dim=c(10,1))
avg8=array(NA,dim=c(10,1))
avg9=array(NA,dim=c(10,1))
avg10=array(NA,dim=c(10,1))
avg11=array(NA,dim=c(10,1))
avg12=array(NA,dim=c(10,1))
avg13=array(NA,dim=c(10,1))
avg14=array(NA,dim=c(10,1))
avg15=array(NA,dim=c(10,1))

sigma1=array(NA,dim=c(10,1))
sigma2=array(NA,dim=c(10,1))
sigma3=array(NA,dim=c(10,1))
sigma4=array(NA,dim=c(10,1))
sigma5=array(NA,dim=c(10,1))
sigma6=array(NA,dim=c(10,1))
sigma7=array(NA,dim=c(10,1))
sigma8=array(NA,dim=c(10,1))
sigma9=array(NA,dim=c(10,1))
sigma10=array(NA,dim=c(10,1))
sigma11=array(NA,dim=c(10,1))
sigma12=array(NA,dim=c(10,1))
sigma13=array(NA,dim=c(10,1))
sigma14=array(NA,dim=c(10,1))
sigma15=array(NA,dim=c(10,1))

for(i in 1:10){
  avg1[i,1]=mean(averagegroup1[i,])
  avg2[i,1]=mean(averagegroup2[i,])
  avg3[i,1]=mean(averagegroup3[i,])
  avg4[i,1]=mean(averagegroup4[i,])
  avg5[i,1]=mean(averagegroup5[i,])
  avg6[i,1]=mean(averagegroup6[i,])
  avg7[i,1]=mean(averagegroup7[i,])
  avg8[i,1]=mean(averagegroup8[i,])
  avg9[i,1]=mean(averagegroup9[i,])
  avg10[i,1]=mean(averagegroup10[i,])
  avg11[i,1]=mean(averagegroup11[i,])
  avg12[i,1]=mean(averagegroup12[i,])
  avg13[i,1]=mean(averagegroup13[i,])
  avg14[i,1]=mean(averagegroup14[i,])
  avg15[i,1]=mean(averagegroup15[i,])
  
  sigma1[i,]=sd(averagegroup1[i,])
  sigma2[i,]=sd(averagegroup2[i,])
  sigma3[i,]=sd(averagegroup3[i,])
  sigma4[i,]=sd(averagegroup4[i,])
  sigma5[i,]=sd(averagegroup5[i,])
  sigma6[i,]=sd(averagegroup6[i,])
  sigma7[i,]=sd(averagegroup7[i,])
  sigma8[i,]=sd(averagegroup8[i,])
  sigma9[i,]=sd(averagegroup9[i,])
  sigma10[i,]=sd(averagegroup10[i,])
  sigma11[i,]=sd(averagegroup11[i,])
  sigma12[i,]=sd(averagegroup12[i,])
  sigma13[i,]=sd(averagegroup13[i,])
  sigma14[i,]=sd(averagegroup14[i,])
  sigma15[i,]=sd(averagegroup15[i,])
}

###### Figure 4: Agent-Based Model Simulation Results
# Only 1st period actual contribution data is used
period<-(1:10)

matplot(avg1[1:10,], type = c("o"),pch=1,col = 1:10,family="Times New Roman",font=1,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 1 over Time: 0.5 multiplier",ann="FALSE",cex.lab=2, cex.axis=2, cex.main=2, cex.sub=2) #plot
polygon(c(period, rev(period)), c(avg1+2*sigma1, rev(avg1-2*sigma1)),
        col = "lightgrey", border = NA)
matlines(trueavg[1,],col="black", lwd=4)
matlines(avg1[1:10,], type = c("o"),pch=1,col = 1:10,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 1 over Time: 0.5 multiplier",ann="FALSE") #plot


matplot(avg2[1:10,], type = c("l"),pch=1,col = 1:10,family="Times New Roman",font=1,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 2 over Time: 0.5 multiplier",ann="FALSE",cex.lab=2, cex.axis=2, cex.main=2, cex.sub=2) #plot
polygon(c(period, rev(period)), c(avg2+2*sigma2, rev(avg2-2*sigma2)),
        col = "lightgrey", border = NA)
matlines(trueavg[2,],col="black", lwd=4)
matlines(avg2[1:10,], type = c("o"),pch=1,col = 1:10,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 2 over Time: 0.5 multiplier",ann="FALSE") #plot


matplot(avg3[1:10,], type = c("o"),pch=1,col = 1:10,family="Times New Roman",font=1,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 3 over Time: 0.5 multiplier",ann="FALSE",cex.lab=2, cex.axis=2, cex.main=2, cex.sub=2) #plot
polygon(c(period, rev(period)), c(avg3+2*sigma3, rev(avg3-2*sigma3)),
        col = "lightgrey", border = NA)
matlines(trueavg[3,],col="black", lwd=4)
matlines(avg3[1:10,], type = c("o"),pch=1,col = 1:10,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 3 over Time: 0.5 multiplier",ann="FALSE") #plot


matplot(avg4[1:10,], type = c("o"),pch=1,col = 1:10,family="Times New Roman",font=1,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 4 over Time: 0.5 multiplier",ann="FALSE",cex.lab=2, cex.axis=2, cex.main=2, cex.sub=2) #plot
polygon(c(period, rev(period)), c(avg4+2*sigma4, rev(avg4-2*sigma4)),
        col = "lightgrey", border = NA)
matlines(trueavg[4,],col="black", lwd=4)
matlines(avg4[1:10,], type = c("o"),pch=1,col = 1:10,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 4 over Time: 0.5 multiplier",ann="FALSE") #plot


matplot(avg5[1:10,], type = c("o"),pch=1,col = 1:10,family="Times New Roman",font=1,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 5 over Time: 0.5 multiplier",ann="FALSE",cex.lab=2, cex.axis=2, cex.main=2, cex.sub=2) #plot
polygon(c(period, rev(period)), c(avg5+2*sigma5, rev(avg5-2*sigma5)),
        col = "lightgrey", border = NA)
matlines(trueavg[5,],col="black", lwd=4)
matlines(avg5[1:10,], type = c("o"),pch=1,col = 1:10,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 5 over Time: 0.5 multiplier",ann="FALSE") #plot


matplot(avg6[1:10,], type = c("o"),pch=1,col = 1:10,family="Times New Roman",font=1,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 6 over Time: 0.5 multiplier",ann="FALSE",cex.lab=2, cex.axis=2, cex.main=2, cex.sub=2) #plot
polygon(c(period, rev(period)), c(avg6+2*sigma6, rev(avg6-2*sigma6)),
        col = "lightgrey", border = NA)
matlines(trueavg[6,],col="black", lwd=4)
matlines(avg6[1:10,], type = c("o"),pch=1,col = 1:10,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 6 over Time: 0.5 multiplier",ann="FALSE") #plot


matplot(avg7[1:10,], type = c("o"),pch=1,col = 1:10,family="Times New Roman",font=1,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 7 over Time: 0.5 multiplier",ann="FALSE",cex.lab=2, cex.axis=2, cex.main=2, cex.sub=2) #plot
polygon(c(period, rev(period)), c(avg7+2*sigma7, rev(avg7-2*sigma7)),
        col = "lightgrey", border = NA)
matlines(trueavg[7,],col="black", lwd=4)
matlines(avg7[1:10,], type = c("o"),pch=1,col = 1:10,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 7 over Time: 0.5 multiplier",ann="FALSE") #plot


matplot(avg8[1:10,], type = c("o"),pch=1,col = 1:10,family="Times New Roman",font=1,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 8 over Time: 0.5 multiplier",ann="FALSE",cex.lab=2, cex.axis=2, cex.main=2, cex.sub=2) #plot
polygon(c(period, rev(period)), c(avg8+2*sigma8, rev(avg8-2*sigma8)),
        col = "lightgrey", border = NA)
matlines(trueavg[8,],col="black", lwd=4)
matlines(avg8[1:10,], type = c("o"),pch=1,col = 1:10,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 8 over Time: 0.5 multiplier",ann="FALSE") #plot


matplot(avg9[1:10,], type = c("o"),pch=1,col = 1:10,family="Times New Roman",font=1,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 9 over Time: 0.5 multiplier",ann="FALSE",cex.lab=2, cex.axis=2, cex.main=2, cex.sub=2) #plot
polygon(c(period, rev(period)), c(avg9+2*sigma9, rev(avg9-2*sigma9)),
        col = "lightgrey", border = NA)
matlines(trueavg[9,],col="black", lwd=4)
matlines(avg9[1:10,], type = c("o"),pch=1,col = 1:10,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 9 over Time: 0.5 multiplier",ann="FALSE") #plot


matplot(avg10[1:10,], type = c("o"),pch=1,col = 1:10,family="Times New Roman",font=1,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 10 over Time: 0.5 multiplier",ann="FALSE",cex.lab=2, cex.axis=2, cex.main=2, cex.sub=2) #plot
polygon(c(period, rev(period)), c(avg10+2*sigma10, rev(avg10-2*sigma10)),
        col = "lightgrey", border = NA)
matlines(trueavg[10,],col="black", lwd=4)
matlines(avg10[1:10,], type = c("o"),pch=1,col = 1:10,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 10 over Time: 0.5 multiplier",ann="FALSE") #plot


matplot(avg11[1:10,], type = c("o"),pch=1,col = 1:10,family="Times New Roman",font=1,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 11 over Time: 0.5 multiplier",ann="FALSE",cex.lab=2, cex.axis=2, cex.main=2, cex.sub=2) #plot
polygon(c(period, rev(period)), c(avg11+2*sigma11, rev(avg11-2*sigma11)),
        col = "lightgrey", border = NA)
matlines(trueavg[11,],col="black", lwd=4)
matlines(avg11[1:10,], type = c("o"),pch=1,col = 1:10,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 11 over Time: 0.5 multiplier",ann="FALSE") #plot


matplot(avg12[1:10,], type = c("o"),pch=1,col = 1:10,family="Times New Roman",font=1,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 12 over Time: 0.5 multiplier",ann="FALSE",cex.lab=2, cex.axis=2, cex.main=2, cex.sub=2) #plot
polygon(c(period, rev(period)), c(avg12+2*sigma12, rev(avg12-2*sigma12)),
        col = "lightgrey", border = NA)
matlines(trueavg[12,],col="black", lwd=4)
matlines(avg12[1:10,], type = c("o"),pch=1,col = 1:10,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 12 over Time: 0.5 multiplier",ann="FALSE") #plot

matplot(avg13[1:10,], type = c("o"),pch=1,col = 1:10,family="Times New Roman",font=1,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 13 over Time: 0.5 multiplier",ann="FALSE",cex.lab=2, cex.axis=2, cex.main=2, cex.sub=2) #plot
polygon(c(period, rev(period)), c(avg13+2*sigma13, rev(avg13-2*sigma13)),
        col = "lightgrey", border = NA)
matlines(trueavg[13,],col="black", lwd=4)
matlines(avg13[1:10,], type = c("o"),pch=1,col = 1:10,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 13 over Time: 0.5 multiplier",ann="FALSE") #plot

matplot(avg14[1:10,], type = c("o"),pch=1,col = 1:10,family="Times New Roman",font=1,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 14 over Time: 0.5 multiplier",ann="FALSE",cex.lab=2, cex.axis=2, cex.main=2, cex.sub=2) #plot
polygon(c(period, rev(period)), c(avg14+2*sigma14, rev(avg14-2*sigma14)),
        col = "lightgrey", border = NA)
matlines(trueavg[14,],col="black", lwd=4)
matlines(avg14[1:10,], type = c("o"),pch=1,col = 1:10,xlab = "Period", ylab = "Average Contribution",ylim=c(-1,21),main="Average Contribution in Group 14 over Time: 0.5 multiplier",ann="FALSE") #plot


matplot(period,avg15[1:10,], type = c("o"),pch=1,col = 1:10,family="Times New Roman",font=1,xlab = "Period", ylab = "Average Contribution", ylim=c(-1,21), main="Average Contribution in Group 15 over Time: 0.5 multiplier", ann="FALSE",cex.lab=2, cex.axis=2, cex.main=2, cex.sub=2) #plot
polygon(c(period, rev(period)), c(avg15+2*sigma15, rev(avg15-2*sigma15)),
        col = "lightgrey", border = NA)
matlines(period,trueavg[15,],col="black", lwd=4)
matlines(period,avg15[1:10,], type = c("o"),pch=1,col = 1:10,xlab = "Period", ylab = "Average Contribution", ylim=c(-1,21), main="Average Contribution in Group 15 over Time: 0.5 multiplier", ann="FALSE") #plot

